function f_min=example_flight(theta,flight)
n=6;% 飞机架数
% theta=zeros(1,n);% 每架飞机调整角初始化
dt=1e-4;% 时间微元
% 本节首先考虑对约束条件不碰撞和区域边缘限制进行软化，通过设计处罚，以使优化
% 结果满足条件限制。
tm=sqrt(2*160^2)/800;% 最大运行时间
v=800;% 速度
flight1=flight;
flight1(:,4)=flight1(:,4)/360*2*pi;% 将角度转换为弧度制
d_crash=zeros(1,n);% 区域内碰撞处罚参数
d_crashm=zeros(2,n);% 飞机到达区域边缘时距离限制处罚参数
for i=0:dt:tm
    for i1=1:n % 模拟飞行过程
        if(flight1(i1,2)>=0&&flight1(i1,2)<=160&&flight1(i1,3)>=0&&flight1(i1,3)<=160)
            flight1(i1,2)=flight1(i1,2)+v*dt*cos(flight1(i1,4)+theta(i1));
            flight1(i1,3)=flight1(i1,3)+v*dt*sin(flight1(i1,4)+theta(i1));
        end
    end
    for i1=1:n % 计算区域内飞机间距离
        if(flight1(i1,2)>0&&flight1(i1,2)<160&&flight1(i1,3)>0&&flight1(i1,3)<160)
            d1=0;
            for i2=1:n
                if(i1~=i2)
                    if(flight1(i2,2)>0&&flight1(i2,2)<160&&flight1(i2,3)>0&& ...
                            flight1(i2,3)<160)
                        d0=sqrt((flight1(i1,2)-flight1(i2,2))^2+(flight1(i1,3)- ...
                            flight1(i2,3))^2);
                        if(d0<8)
                            d1=d1+1/(d0+eps)*dt/tm;
                        end
                    end
                end
            end
            d_crash(i1)=d_crash(i1)+d1;
        end
    end
    for i1=1:n % 计算飞机在区域边缘时情况
        if(flight1(i1,2)<=0||flight1(i1,2)>=160||flight1(i1,3)<=0||flight1(i1,3)>=160)
            if(d_crashm(1,i1)==0)
                d1=0;
                for i2=1:n
                    if(i1~=i2)
                        if(flight1(i2,2)>0&&flight1(i2,2)<160&&flight1(i2,3)>0&& ...
                                flight1(i2,3)<160)
                            d0=sqrt((flight1(i1,2)-flight1(i2,2))^2+(flight1(i1,3)- ...
                                flight1(i2,3))^2);
                            if(d0<60)
                                d1=d1+60/(d0+eps);
                            end
                        end
                    end
                end
                d_crashm(2,i1)=d_crashm(2,i1)+d1;
                d_crashm(1,i1)=1;
            end
        end
    end
end
rho1=100;rho2=100;% 处罚因子
f_min=sum(abs(theta))/pi*180+rho1*sum(d_crash)+rho2*sum(d_crashm(2,:));% 目标函数设计